Benchmarking Defmod, an open source FEM code for modeling episodic fault rupture

We present Defmod, an open source (linear) finite element code that enables us to efficiently model the crustal deformation due to (quasi-)static and dynamic loadings, poroelastic flow, viscoelastic flow and frictional fault slip. Ali (2015) provides the original code introducing an implicit solver for (quasi-)static problem, and an explicit solver for dynamic problem. The fault constraint is implemented via Lagrange Multiplier. Meng (2015) combines these two solvers into a hybrid solver that uses failure criteria and friction laws to adaptively switch between the (quasi-)static state and dynamic state. The code is capable of modeling episodic fault rupture driven by quasi-static loadings, e.g. due to reservoir fluid withdraw or injection. Here, we focus on benchmarking the Defmod results against some establish results. 1. Quasi-static crustal deformation When a region is subjected to a gradual loading process, such as tectonic stress changes, viscoelastic relaxation, and pore pressure changes, it deforms in a quasi-static manner. Every snapshot of a quasi-static process, as opposed to a dynamic process, satisfies stress equilibrium. The inertial force is considered negligible, since the net force is small, and the time scale is large. For linear constitutive law and small strain problems, the finite element method, Zienkiewicz (2000), provides a system of linear equations describing the (quasi-)static state. Eq. (1) lists the absolute and incremental versions of the linear equation, Smith and Griffiths, 2004. Δ K U F K U F = , absolute, = Δ , incremental. n n n n n n (1) where, K is the system stiffness matrix, U is the solution vector and F is the nodal force, including fluid source. The subscript n is the time index. In this study, we use the incremental equation. The solution space U Δ n of a poroelastic problem contains the nodal displacement and pressure, ⎡ ⎣⎢ ⎤ ⎦⎥ U Δ = n u p Δ Δ n n , whereas the solution for an elastic problem is only the nodal displacement. The stiffness matrix Kn and RHS function F Δ n are also different for the elastic and poroelastic problems, Eq. (2). ⎧ ⎨⎪ ⎩⎪ ⎛ ⎝ ⎜⎜ ⎡ ⎣⎢ ⎤ ⎦⎥ ⎡ ⎣⎢ ⎤ ⎦⎥ ⎞ ⎠ ⎟⎟ t t K F K f K H H K S f q K p ( , Δ ) = ( , Δ ) , elastic. − Δ + , Δ − Δ , poroelastic. n n e n e T c p n n c n−1 (2) where, Ke is the elastic stiffness matrix, depending on the elastic constants of the solid. Kc is the fluid stiffness matrix, depending on the fluid flow conductivity. H is the coupling matrix, depending on the Biot's coefficient. Sp is the storage matrix, depending the solid compressibility and porosity, and the fluid compressibility. Smith and Griffiths (2004) provide the detailed formulation for these matrices and vectors. Note, the stiffness matrix Kn is constant for evenly spaced time step t Δ . In a later section, we show that for Newtonian viscoelasticity, Kn is, although modified, still independent of time. fn and qn are nodal force and fluid source respectively. The detailed formulations of these matrices and RHS vectors are given in Appendix A 2. Poroelastic model and benchmark Unstable pressure is caused by using linear elements, known as the Ladyzenskaja-Babuska-Brezzi restrictions. The local pressure projection scheme, Bochev and Dohrmann (2006), is implemented to stabilize the pore pressure, ⎡ ⎣⎢ ⎤ ⎦⎥ ⎡ ⎣⎢ ⎤ ⎦⎥ ∫ n I n I G dΩ F F 0 H p K K 0 0 0 H H N N = − , = + , where , = ( − (1/ ) ) ( − (1/ ) )/(2 ) , n n s n n n s s Ω e e N N +1 +1 +1 +1